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Abstract 

A physical model and interactive computer graphics techniques have been developed for the vi- 
sualization of the basic physical process of stochastic dispersion and mixing from steady-state 
CFD calculations. The mixing of massless particles and inertial particles is visualized by trans- 
forming the vector field from a traditionally Eulerian reference frame into a Lagrangian reference 
frame. Groups of particles are traced through the vector field for the mean path as well as their 
statistical dispersion about the mean position by using added scalar information about the root 
mean square value of the vector field and its Lagrangian time scale. In this way, clouds of par- 
ticles in a turbulent environment are traced, not just mean paths. In combustion simulations of 
many industrial processes, good mixing is required to achieve a sufficient degree of combustion 
efficiency. The ability to visualize this multiphase mixing can not only help identify poor mixing 
but also explain the mechanism for poor mixing. The information gained from the visualization 
can be used to improve the overall combustion efficiency in utility boilers or propulsion devices. 
We have used this technique to visualize steady-state simulations of the combustion performance 
in several furnace designs. 


iThis research was supported in part by the National Aeronautics and Space Administration under NASA 
contract NAS1-19480 while the author wets in residence at the Institute for Computer Applications in Science and 
Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23681-0001. 




1 Introduction 


In computational studies of reactive-flow applications, such as combustion and chemical reac- 
tor design and trouble-shooting, the physical mechanisms of dispersion and mixing are often 
what the engineer is trying to understand. For the results of these computations, traditional 
visualization methods such as particle tracing, cutting planes and contour surface plots do not 
provide this mechanistic insight since they only display the result of mixing, they do not display 
the process of mixing. A cloud tracing technique has been developed to address visualization of 
these mechanisms [8]. In this paper, we extend that method to include inertial particles and to 
clarify the physical basis and properties of the visualization of these clouds. 

The cloud visualization method is based on transforming the vector field from a traditionally 
Eulerian (fixed) reference frame into a Lagrangian (moving) reference frame. Mean fluid elements 
are traced through the computed vector field (as in traditional particle tracing methods). The 
statistical dispersion of the fluid elements about the mean position is also computed by using 
added scalar information about the root mean square value for the vector field and its Lagrangian 
time scale. These scalar values can usually be computed easily from existing information in 
computational fluid dynamics codes (i.e. turbulence models, etc.). In this way ensemble averaged 
“clouds of particles” are traced not just mean paths. Clouds of fluid elements and inertial particles 
are useful in distinguishing multi-phase mixing effects. 

Each cloud follows the mean particle path and has a tube like shape, which should be distin- 
guished from the traditional streamtubes [2] and their variants [5]. In essence, clouds differ from 
streamtubes in that they are derived from mean components of the velocity vector but come 
from fluctuating components. A cloud can cross streamlines due to eddy transport whereas a 
streamtube does not. A streamtube may be treated as if isolated from the adjacent fluid. As 
will be shown, a cloud can be viewed as the ensemble average, spatial probability distribution 
function of all particles starting from a given location. 

The cloud visualization procedure allows a user to independently select what scalar value 
should be used to represent the color of the surface of a cloud. This is particularly useful in 
multi-phase mixing and reaction applications since the scalar properties of either phase can 
be used to color the cloud. The properties of the cloud are computed from ensemble average 
statistics of the dispersing particles due to fluid turbulence. This characteristic is discussed by 
Jain [7] and results in a visualization algorithm that is closer to the physical mixing process than 
alternate visualization schemes. 

The physical model used to simulate the turbulent dispersion phenomenon is described in 
Section 2. Visualization results derived based the model are shown and discussed in Section 3. 
The three-dimensional computer graphics techniques that we have developed for this research 
are described in [8]. 
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2 A 3D Turbulent Particle Dispersion Model 

The dispersion of particles or fluid elements by turbulent motion has been an active field of 
research in recent years. Practical interests are widespread and include the combustion of partic- 
ulate and droplet fuels, the pneumatic transport of particles, the mixing and settling of particles, 
particle deposition, electrostatic precipitation, and particle fluid separations. Practical com- 
bustion processes like spray combustion, pulverized coal and coal slurry combustion, fluidized 
beds, sorbent injection, and hazardous waste incineration, introduce the fuel as particles, liquid 
droplets, or slurries into turbulent environments. 

Mixing and dispersion of gaseous streams involves fluid elements that disperse like particles but 
without any added inertial effects of particles with drag. The equations used here were developed 
by Batchelor [1]. They have been extended and used for turbulent dispersion of inertial particles 
in a model for the stochastic transport of particles, referred to as the STP model [3]. This 
approach is an efficient stochastic Lagrangian scheme that does not compromise accuracy. A 3D 
implementation and extension of the scheme has been described by Jain [7]. 

Mathematically, in the particle cloud approach, we are tracking the evolution of a probability 
density function (pdf) for particle position with respect to residence time. A multivariate normal 
pdf is used to describe the distribution of particles in the cloud. The general formula is 


P(xi, t) 
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However, there is rarely enough information about a turbulent flow field to evaluate the off 
diagonal terms of the variance tensor cr,j. If they are assumed to be zero, the resulting pdf has 
the following form: 


P(Xi, t ) 


1 

(87t)( 3 / 2 ) nf=i 



( 3 ) 


where 
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In the above equations the pdf is completely described in terms of two vectors, the mean (expected 
value), Hi, and the variance (second moment), a 2 *, which vary explicitly with residence time. 

*The notation of for the variance in the i direction refers to the ca component of the variance tensor. 
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2.1 Mean of the Particle pdf 

The mean of the pdf at some residence time, represents the most likely location of the particles 
in the cloud. The mean location is obtained by integrating the particle velocity which is defined 
by an equation of motion for the cloud particles. 

m(t) = ( xi(t )) = [ (K(*i))<fti + (*,-(0)) (5) 

JO 

For a small heavy spherical particle the Basset-Bousineq-Oseen (BBO) equation of motion re- 
duces to the form [9]. 

m ’^dT = "’H&W " (K)) + (m " - ro /)»< (6) 

The first term on the right hand side of the equation represents the force acting on a particle 
due to the aerodynamic drag. The drag coefficient, is calculated using a non-linear Stokes law, 
in terms of the particle Reynolds number, 

p fi f 

In the creeping flow regime the drag coefficient can be obtained by using a non-linear Stoke’s 
law relationship 

CD = (^) (1+0 - 15fle °' 6 ‘ 7) (*) 

The second term on the right hand side of Eq. 6 denotes external forces such as gravity. Additional 
terms in the equation of motion, describing virtual mass, thrust, Magnus, Basset, and Saffman 
forces, etc. have a contribution of less than 1% for the particles under consideration in various 
flows. Consequently, they are neglected in this analysis. 

The angle brackets, {), used in the equation of motion imply ensemble average values. If N 
identical experiments are performed, then an ensemble average is defined as, 

W = (9) 

The ensemble average represents the most likely occurrence of an event with several realizations. 
Thus, the particle trajectory obtained, by solving an equation of motion in terms of ensemble 
average values, has the weighted effect of all the possible realizations of the particle(s) represented 
by the cloud and represents the statistical mean of the cloud. The ensemble average value can 
be equated mathematically to the expected value or the first moment of a positional pdf, which 
is the basis for the first equality in Eq. 5. An ensemble average gas velocity can be obtained as 
the expected value of its pdf as, 

(U(xi,t)) = r U(x,t)P(U)dU (10) 

J — oo 
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2.2 Variance of the Particle pdf 

The other variable required for defining the cloud is the variance, a 2 which is a measure of 
the spread of the pdf around the mean location in each co-ordinate direction. The variance is 
defined in terms of the fluctuating components. If we write the instantaneous values of velocity 
and position in terms of the mean and fluctuating components as, 

w = (V) + t; / (11) 


x=(X)+ x' 


( 12 ) 


then the variance of the pdf, by definition, is the ensemble average of the mean square displace- 
ments of the cloud particle positions. The ordinary differential equation governing its evolution 
with time can be derived as: 

^ ~ = 2 J‘mw*))dh (13) 

This result can be expressed in terms of a classical turbulence property, namely, the particle 
velocity correlation function. It’s basis lies in Taylor’s work on the Lagrangian description of 
particle dispersion and is defined as [6] 






(14) 


Substituting this in Eq. 13 and integrating provides a relationship for the variance in terms 
of two particle turbulence statistical quantities: the mean square velocity fluctuations, and the 
particle velocity correlation function, R p (t ), as originally derived by Batchelor [1]: 


aij(t) = 2 [ [ (u-(ti)u'(t 2 ))dtidt 2 

Jo Jo 


(15) 


<r«(i) = 2 /' /‘ ! [(ff((.)<(( 2 ))]’ /2 K(l 7,h) + RUh, (16) 

J 0 J 0 

The variance tensor, cr,j, thus obtained, would be used in the general form of the pdf given by 
Eq. 1 and Eq. 2. 

2.3 The Particle Velocity Correlation Function 

The particle velocity correlation as defined by Eq. 14 represents the rate at which the instanta- 
neous velocity of the particle at some point (corresponding to residence time, t) loses its effect as 
we move away (to a location corresponding to time <i). It is a measure of the time scale, T, over 
which a velocity fluctuation remains correlated. At r = 0 the velocity is fully correlated with 
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itself and its value is 1. As we move away, the particle velocity loses its correlation monotoni- 
cally and drops asymptotically to 0. The particle velocity correlation, R p , is different from the 
Lagrangian fluid velocity correlation, R } . which can be defined similarly, in terms of the fluid 
fluctuating velocities. 

The equations presented in this section have been derived by Wang [12] and further mathe- 
matical details of the derivation can be found in his dissertation. 

The derivation proceeds in two main steps. First, the particle velocity correlation is related 
to that of a fluid element moving along the particle trajectory. The correlation function for 
this fluid element is then related to fluid Eulerian velocity correlations. Combining the above 
two relations the particle velocity correlation is found in terms of three dimensionless Eulerian 
turbulence parameters; the Stokes number, the drift parameter, and the turbulence structure 
parameter. Since the fluctuating components of the velocities are of interest here, this analysis 
is developed in a moving Eulerian reference frame. 

Most analyses in the past have not distinguished between the particle and fluid element cor- 
relations. Reeks [11], and Pismen and Nir [10] have taken this into account. However, their 
complex analysis of the relation between fluid Lagrangian and Eulerian velocities has made their 
final results mathematically unwieldy. Wang has avoided this by using an approach similar to 
Csanady’s [4], while including inertial effects at the same time. 

The particle velocity correlation tensor 


RU t ) = WWi(t + r)> 

is related to the fluid velocity correlation tensor along the particle trajectory, namely 
R !j( T ) = <«-(0)uJ(t + r) = (u[(x(0),0)u'(x(r),r))) 

The particle equation of motion can be integrated [11] to give the particle velocity 


f* U.-(x(t),t) ,T-i N , c 
exp( )dr -f v d S i3 


(17) 


(18) 


(19) 


Maxey [9] has shown that the mean velocity of a particle is slightly larger than v d (by less 
than 3-5% in most cases). If this difference is neglected the instantaneous fluctuating velocity of 
the particle can be found as 


= Vi(t) - ( Vi(t )) = Vi(t) - V d Si 3 

Combining Eq. 18, Eq. 19, and Eq. 20 

RP ij( T ) = -K / / R ij( T i - T 2 )exp( — +T 2 — -)dr 0 dr 2 

~ a J- oo J — oo T a 

Partial integration of Eq. 21 gives 

R U t ) =97 / ex P( — ) R ij U)d£ 

£T a J — oo T a 


( 20 ) 


( 21 ) 


( 22 ) 
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This relation shows the difference between the velocity correlation for a particle and a fluid 
element moving along the particle trajectory, due to particle inertia. Reeks [11] has derived this 
same equation. It can be seen that Csanady’s assumption of RP = R* is good only if r a — >0, i.e., 
zero inertia particles. 

Next, to relate R?(t) to available fluid turbulence statistics, Csanady’s approach [4] is used 
and modified to include the effects of inertia in the time scale. In this approach, the normalized 
velocity autocorrelations of the turbulent flow in the moving Eulerian frame are assumed to be 
exponential functions, namely, 

D l {t) = exp(-^r) (23) 


D mE (T) = exp(--P-) 
J-mE 


f(r) = exp(-y-) 
L S 


(24) 

(25) 


9(f) = (1 - ^)e*p(- jj) ( 26 ) 

Dl is the Lagrangian fluid velocity correlation and D mE is the one-point Eulerian fluid velocity 
temporal correlation. f(r ) and g(r) are the longitudinal and transverse fluid velocity spatial 
correlations, respectively. Exponential correlations have been shown to fit experimental data [6]. 

For the fluid velocity correlation along the particle trajectory, with time delay, t, the mean 
spatial displacement of the particle in this time would be Vdt. Hence, a spatial correlation should 
be included in R{j(r) to account for the relative displacement. In the limiting case of v<i—> oo, 
the particle displacement is relatively large and the spatial correlation is dominant. Thus, we 
have R{ x = /?{ 2 (r) = g(vdT ) and R1 3 (t) = /(u^r). In the other extreme, when u<£— >0, i2-(r) is 
mainly a temporal correlation. The time scale T, for this correlation is a function of the particle 
inertia, and is not the same as Ti. 

In general, when the effects of both inertia and drift velocity are present, there should be a 
smooth transition of R{j(r) from a fluid temporal velocity correlation to a fluid spatial velocity 
correlation, when the drift velocity increases from 0 to oo. This is done by using Csanady’s 
hypothesis that R& (r) is constant on the ellipses 

r 2 V?T 2 

— + -p- = constant ( 27 ) 

Using the velocity correlations Eq. 23 - Eq. 27 , R{j(T) is obtained as 

-RiiM = RL(t) = « 2 ( 1 ~ |^)exp(-^\/l + mr 7 2 ) ( 28 ) 


6 



TV , x 

m-T = — — (31) 

Substituting Eq. 28 or Eq. 29 into Eq. 22, the following relationship for the particle velocity 
correlation is obtained 


R h ( T ) = {-Q-) ex P{--^r){StT\] 1 + " 4 7 2 - + ^ r7 - + 1 ) + 



+ mTSth ^ + + o.5^) 

(32) 


W = ^tr^ff exp{ _ r } _ 

(33) 

where, St is 

the Stokes number and is a measure of the particles inertia. 



II 

tx 

(34) 

and 

0 = Stj( 1 + m\ 7 2 ) — 1 

(35) 


The time scale T, which is physically the time scale for a fluid element traveling along the 
particle trajectory, is needed. It is a function of the particle inertia. In the limit that St~^> 1, 
the particle’s inertia is very large, so that the particle responds very slowly to fluctuations in 
fluid velocity. Consequently, the particle stays at one point for a time interval of the order of 
T m E . Thus when St^> 1. On the other hand if ST— >0, the particle reduces to a fluid element and 
T— >Tl. In general T is a function of St and varies between Ti and T mE . Wang and Stock [13] 
have done a numerical simulation using turbulence generated by Fourier modes and obtained the 
following curve fit for T as a function of the Stokes parameter. 

T(St) _ i constant /nr%x 

T mE = ” (1 + ST)°- 4 ( 1+0 - 01St ) ( 36 ) 

where, the constant=l - {Ti/TmE) = 0.644, for m = 1. 
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To get the dimensionless parameters Stx and mj, the parameters St and m based on Eulerian 
time scales are first calculated as, 


St = 


and 


m = 


T m E 

TmEU 

Lj 


Eq. 36 gives T(St). St? and mx are then calculated as 

T m E 


St T = St 


and 


mx = m 


T(St) 

TjSt ) 
T m E 


We calculate T m E from the fixed Eulerian parameters by an approximate relation as 

T m E = TfE~ 
u' 

Tl is related to T m E by the turbulence parameter, m. 

Tl = G{rn)T m E 


(37) 

(38) 

(39) 

(40) 

(41) 

(42) 



exp(-y 2 ) 

1 + f[V*y er f(y) - 1 + exp(-y 2 )] 


dy 


(43) 


G(m) goes from 1 to 0 asymptotically as m increases. This trend follows our discussion earlier, 
wherein as the turbulence increases, the fluid element is likely to escape the eddy quicker and Tl 
will decrease. 

The behavior of the correlation function for various Stokes and drift parameters has been 
studied [12] and shown to realistically describe the different effects of particle behavior. The 
correlations show a smooth transition between the drift velocity and inertial effects. Also, these 
correlations are anisotropic providing greater dispersion in the direction parallel to the drift 
velocity as per the continuity effect. This anisotropy can occur even if the fluid phase turbulence 
statistics available are isotropic (from the k — e turbulence model, for example). This effect can 
be important in determining the cloud position in complex mixing systems. Wang has shown 
parametric zones where particle dispersion is more than the fluid dispersion as has been reported 
in the past. These correlations incorporate a physically realistic representation of the inertial 
particle effects. 
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3 Visualization Results 


We have used the design of a low emission pulverized-coal combustion boiler as a demonstrative 
problem to show the utility of this visualization method for mixing of gases and particles. In 
this slagging boiler design the pulverized-coal fuel is feed through 10 separate burners into a 
U-shaped furnace. Figure la shows a rendered image of the furnace geometry revealing the U- 
shaped structure of the furnace. This shape allows for the ash to be collected as a flowing slag 
and ’’tapped” at the bottom of the U. The air is staged into three separate levels of injection 
to help in minimizing pollutant formation. This design is being considered for the U.S. DOE 
Combustion 2000 program for the design of a low emissions boiler system (LEBS) for the next 
century. The combustion computations have been performed by one of the authors (Smith) in 
support of this DOE project. The design question is to determine from numerical simulations 
whether the air sufficiently mixes with the pulverized-coal stream or not. We have used the cloud 
tracing visualization of both the fluid (air) streams and the inertial (pulverized-coal) streams to 
give insight into the extent of mixing. 

The images shown in Figures lb and lc display some aspects of flow and reaction within 
the boiler. These are direct volume rendered images of the CO concentration (lb) and 0 2 
concentration (lc) in three-dimensional space. The view is from the side of the furnace. This more 
traditional visualization tool gives the combustion engineer an idea about important operational 
variables and infers something about the mixing process that occurred to create the results shown 
in this image. Before it reacts, the fuel (as visualized by the CO concentration in Figure lb) 
mixes with the air (as visualized by the 0 2 concentration in Figure lc) to create regions of 
higher temperature. However, we still do not see where the particle streams went, or where the 
air streams went, or how the different streams interact. 

Figures Id through lg use the cloud dispersion visualization scheme of this paper. Figures Id, 
le and If compare the mixing and dispersion of the entraining air of part of one of the pulverized- 
coal streams (If) with the coal particles from the same part of the same stream (Id and le). Coal 
particle clouds of two different particle sizes (150mm in Id and 10mm in le) are compared with 
the fluid dispersion (If). The surface of the dispersing cloud is colored by two different scalars 
(temperature for Id &: le and 0 2 concentration for If) in these images. The fluid (air) streams 
which disperse most rapidly is obvious from this visualization. The smallest particles disperse 
similarly to the gases but more slowly. The largest particles disperse more slowly. The inertial 
particles also interact with various surfaces in the furnace and are partially or totally collected or 
removed from the cloud. Thus the clouds also diminish in size at various locations. One location 
in particular is just downstream from the bottom of the U in this furnace. At this location there 
is a slag screen consisting of a bank of cooled tubes crossing the furnace at an angle across the 
full furnace. Here many of the particles that have not impacted on the slagging furnace walls are 
collected on the screen. Some particles penetrate this screen and continue through the furnace 
section. 

Figures lg through lh use transparency with the rendered image to show the cloud dispersion 
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of the same clouds as shown in Figures If but also includes clouds initiated in air jets starting 
on the side of the furnace, which allows us to see the interaction of these horizontal and vertical 
jets. All the figures to this point have shown the cloud surface at a distance of one standard 
deviation from the mean location of the cloud. Figures lg-lh use 0.34 standard deviations to 
emphasize the differences between the two streams. Figure lh shows a close up view of Figure 
lg. The surface in both images is colored by the scalar value of the O 2 concentration. Figure li 
shows similar views of another set of interacting clouds. 

All of these cloud tracing images show the extent of mixing of air with the reacting coal particle 
stream. This mixing determines the efficiency of the low emission boiler concept. The properties 
of the cloud in this visualization method are determined from the ensemble averaged variables 
across the spatial distribution at each mean residence time of the cloud. 

The cloud surfaces in all of these visualizations are represented as triangular meshes which 
can be efficiently rendered by conventional graphics hardware. The surface is colored by using 
cloud properties selected by the user. During a visualization session, the user can interactively 
pick seed positions in the data domain and see how flows starting from these positions develop 
and mix together in some way. 

4 Conclusions 

We have shown that by using cloud dispersion concepts it is possible to visualize the dispersion 
and mixing of both massless or fluid particles and inertial particle systems. The cloud can inter- 
act with surfaces by sticking or rebounding the particles within the cloud. The use of hardware 
surface rendering provides rapid feedback on mixing and dispersion processes, and is computa- 
tionally efficient enough to be interactive. This visualization technique is an effective analysis tool 
for practical engineering problems by applying it to data sets from three-dimensional combustion 
simulations. Combustion involves many physical and chemical processes whose fundamentals are 
not yet completely understood. This technique provides an improved insight into the mecha- 
nisms controlling high-speed mixing and combustion in utility boilers or propulsion devices. It 
gives visualization to the mechanism of mixing not just the result of mixing. Further insight can 
be obtained by using volume rendering so that the cloud content can also be visualized, not just 
the exterior surface of the cloud. 
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Figure 1 : Visualization of Gas and Particle Mixing 
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5 Nomenclature 


Symbol 

Units 

Definition 

Cd 

__-3 

771 

drag coefficient in the particle equation of motion 

D l 


Lagrangian temporal fluid velocity correlation function 

D m E 


One-point temporal fluid velocity correlation function 

Lf 

m 

typical eddy length scale for fluid turbulence 

m p 

kg 

mass of an individual particle 

m 


turbulence structure parameter 

n 

1 

771 3 

average number density of particles 

P 

1 

m 3 

probability density function 

BP 


particle velocity correlation function 

St 


Stokes number 

T g 


time scale of gas turbulence 

T fE 

5 

correlation time scale in a fixed Eulerian reference frame 

T m E 

S 

correlation time scale in a moving Eulerian reference frame 

T l 

S 

correlation time scale in a Lagrangian reference frame 

T 

s 

particle correlation time 

t 

s 

residence time 


s 

arbitrary time 

U 

771 

s 

instantaneous gas velocity 

u 

m 

s 

average gas velocity 

u' 

m 

s 

gas velocity fluctuation 

Vtd 

m 

s 

turbulent diffusion velocity of the particle 

V 

m 

s 

instantaneous particle velocity 

V 

771 

S 

average particle velocity 

v' 

m 

s 

particle velocity fluctuation 

X 

m 

particle location 

X 

m 

average particle location 

x' 

771 

fluctuation in particle position 

x,y,z 

771 

coordinate length 

0 

771 3 •$ 

momentum exchange coefficient between gas and particle 

7 


drift parameter 

0 


constant defined by Equation 35 

P 

m 

mean of particle position pdf 

Pf 

777*5 

laminar viscosity of fluid phase 

PS 

JSL 

m 3 

density of fluid phase 

Pp 

m 3 

particle density 

a 

771 

standard deviation of pdf for particle position 

Sij 

- 

Kronecker delta 
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